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ABSTRACT 

We present a reduction algorithm that simultaneously ex- 
tends Hermite's reduction for rational functions and the 
Hermite-like reduction for hyperexponential functions. It 
yields a unique additive decomposition and allows to de- 
cide hyperexponential integrability. Based on this reduction 
algorithm, we design a new method to compute minimal 
telescopers for bivariate hyperexponential functions. One of 
its main features is that it can avoid the costly computa- 
tion of certificates. Its implementation outperforms Maple's 
function DEtools[Zeilberger]. Moreover, we derive an order 
bound on minimal telescopers, which is more general and 
tighter than the known one. 

Categories and Subject Descriptors 

1.1.2 [Computing Methodologies]: Symbolic and Alge- 
braic Manipulation — Algebraic Algorithms 

General Terms 

Algorithms, Theory 

Keywords 
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1. INTRODUCTION 

Given a univariate rational function r, Hermite reduc- 
tion in [12, 13, 6] finds rational functions r\ and ri s.t. 
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(i) r — ri+r2, (ii) ri is rational integrable, (iii) r% is a 
proper fraction with a squarefree denominator. The addi- 
tive decomposition is unique, and r is rational integrable if 
and only if r2 = 0. 

A univariate function is hyperexponential if its logarith- 
mic derivative is rational. Exponential, radical and rational 
functions are hyperexponential. Rational Hermite reduction 
has been extended to hyperexponential functions by Daven- 
port in [9] and by Geddes, Le and Li in [10]. The former 
aims at solving Risch's equation; the latter is a differen- 
tial analogue of the reduction algorithm for hypergeometric 
terms in [2]. For a given hyperexponential function H, the 
reduction algorithms in [9, 10] compute two hyperexponen- 
tial functions Hi and H2 s.t. (i) H = Hi + H2, (ii) Hi 
is hyperexponential integrable, (iii) H2 is minimal in some 
sense. However, H2 is not unique in general and it may be 
nonzero even when H is hyperexponential integrable. In or- 
der to decide the integrability of H, one additionally needs 
to compute polynomial solutions of a first-order linear dif- 
ferential equation. 

The method of creative telescoping for hyperexponential 
functions is developed by Almkvist and Zeilberger in [3]. It 
is nowadays an important automatic tool for computing def- 
inite integrals. Recently, it has also played an important role 
in the resolution of intriguing problems in enumerative com- 
binatorics [14, 15]. For a bivariate hyperexponential func- 
tion H(x,y), the problem of creative telescoping is to find a 
nonzero operator L(x, D x )£¥(x){D x ) , the ring of linear dif- 
ferential operators over the rational-function field ¥(x), s.t. 



L(x,D x )(H) = D y (G) 



(1) 



for some hyperexponential function G, where D x = d/dx 
and D y — d/dy. The operator L above is called a telescoper 
for H, and G is the corresponding certificate. An algorithm 
for solving (1) is given in [3], and is based on differential 
Gosper's algorithm. An algorithm for rational-function tele- 
scoping is given in [5], and is based on Hermite reduction. 
The latter separates the computation for telescopers from 
that for certificates, and has a lower complexity than the 
former for rational functions. 

In the present paper, we develop a reduction algorithm 
which, given a univariate hyperexponential function H , con- 
structs two hyperexponential functions Hi and H2 s.t. (i) 
H — Hi + H2, (ii) Hi is hyperexponential integrable, and 



(iii) H2 is either zero or not hyperexponential integrable. We 
show that H2 in the above additive decomposition is unique 
and can be obtained without computing polynomial solu- 
tions of any differential equation. Our algorithm is based 
on the Hermite-like reduction in [10], a differential variant 
of the polynomial reduction in [2] and on the idea for re- 
ducing simple radicals in [18, Proposition 7]. The main new 
ingredients are the uniqueness of H2 and an easy way to 
compute H2, which are crucial for many applications. These 
enable us to extend the reduction-based rational telescoping 
algorithm in [5] to the hyperexponential case, and derive an 
order bound on the telescopers. The bound is more general 
and tighter than that given in [4]. 

The rest of the paper is organized as follows. We review 
the notion of hyperexponential functions and Hermite-like 
reduction in Sections 2 and 3, respectively. A new reduc- 
tion algorithm is developed for hyperexponential functions 
in Section 4. After introducing kernel reduction in Sec- 
tion 5, we present a reduction-based telescoping algorithm 
for bivariate hyperexponential functions, and derive an up- 
per bound on the order of minimal telescopers in Section 6. 
We briefly describe an implementation of the new telescop- 
ing algorithm, and present some experimental results in Sec- 
tion 7, which validate its practical relevance. 

As a matter of notation, we let E be a field of characteristic 
zero and E(y) be the field of rational functions in y over E. 
For a polynomial p £ E[y], we denote by deg(p) and lc(p) 
the degree and leading coefficient of p, respectively. Let D y 
denote the usual derivation d/dy on E(y). Then (E(y),D y ) 
is a differential field. 

2. HYPEREXPONENTIAL FUNCTIONS 

Hyperexponential functions share the common properties 
of rational functions, simple radicals, and exponential func- 
tions. Together with hypergeometric terms, they are fre- 
quently viewed as a special and important class of "closed- 
form" solutions of linear differential and difference equations 
with polynomial coefficients. 

Definition 1. Let $ be a differential field extension o/E(y). 
A nonzero element H £ $ is said to be hyperexponential 
over E(y) if its logarithmic derivative D y (H)/H £ E(y). 

The product of hyperexponential functions is also hyper- 
exponential. Two hyperexponential functions H\,H2 are 
said to be similar if there exists r £ E(y) s.t. H 1 = rH2- 
The sum of similar hyperexponential functions is still hy- 
perexponential, provided that it is nonzero. 

For brevity, we use the notation exp(J fdy) to indicate a 
hyperexponential function whose logarithmic derivative is /. 
For a rational function r £ E, we have 

rexp (/^ dy ) =exp (/ (Z + A/MA') dy) ■ 

A univariate hyperexponential function H is said to be hy- 
perexponential integrable if it is the derivative of another hy- 
perexponential function. For brevity, we say "integrable" in- 
stead of "hyperexponential integrable" in the sequel. 

Assume that H — r exp (f fdy) is integrable. Then H 
is equal to D y (G) for some hyperexponential function G. A 
straightforward calculation shows that G is similar to D y (G), 
and so is H. Set G = sexp(J/) for some s £ E(y). 
Then H = D V (G) if and only if 

r = D y {s) + fs. (2) 



Deciding the integrability of H amounts to finding a rational 
solution s s.t. the above equation holds. 

3. HERMITE-LIKE REDUCTION 

Reduction algorithms have been developed for computing 
additive decompositions of rational functions [16, 12, 13], hy- 
pergeometric terms [1, 2], and hyperexponential functions [9, 
10]. Those algorithms can be viewed as generalizations of 
Gosper's algorithm [11] and its differential analogue [3, §5]. 

For a hyperexponential function H, a reduction algorithm 
computes two hyperexponential functions Hi,H2 s.t. 

H = D y (Hi) + H 2 . (3) 

It turns out that H, Hi and H2 are similar. So we may 
write H=rexp (J /dy) and Hi=Vi exp (J /dy) , where r, n, f 
belong to E(y) and i = 1, 2. Then (3) translates into 

r = D y (n) + fn +r 2 . 

A reduction algorithm for computing (3) amounts to choos- 
ing rational functions r, / and r\ so that r2 satisfies prop- 
erties similar to those obtained in Hermite reduction for ra- 
tional functions. There are at least two approaches to this 
end. One is given in [9], and the other in [10]. We review 
the latter, because the notion of differential-reduced rational 
functions plays a key role in Lemma 6 in Section 4. 

Recall [10, §2] that a rational function r — a/b £ E(y) is 
said to be differential-reduced w.r.t. y if 

gcd (fe, a — i D y (b)) = 1 for all i £ Z. 

By Lemma 2 in [10], r is differential-reduced if and only 
if none of its residues is an integer. The differential ra- 
tional canonical form of a rational function / in E(y) is a 
pair (K, S) in E(y) x E(y) s.t. (i) K is differential-reduced; 

(ii) the denominator of S is coprime with that of K; and 

(iii) / is equal to K + D y (S) / S. Every rational function has 
a unique canonical form in the sense that K is unique and S 
is unique up to a multiplicative constant in E [10, §3]. We 
call K and S the kernel and shell of /, respectively. They 
can be constructed by the method described in [10, §3]. 

Let H be a univariate hyperexponential function in the 
form exp( J fdy) over E. Assume that K and S are the kernel 
and shell of /, respectively. Then H = S exp (J K dy) . Note 
that K = if and only if H is a rational function, which is 
equal to cS for some c £ E. 

Example 2. Let H = \/y 2 + l/(y — l) 2 - The logarithmic 
derivative of H is 

D y H _ D y (l/(y-l) 2 ) y 

H i/(y-i) 2 y 2 + i' 

where y/(y 2 + l) is differential-reduced. The kernel and shell 
of D y (H)/H are y/(y 2 + 1) and l/(y — l) 2 , respectively. 
So H = exp (/ y/(y 2 + 1) dy) /(y - l) 2 . 

For brevity, we make a notational convention. 

Convention 3. Let H denote a hyperexponential function 
whose logarithmic derivative has kernel K and shell S. As- 
sume that K is nonzero, that is, H is not a rational func- 
tion. Set T — exp (J K dy). Moreover, write K = k\/k2, 
where ki, are polynomials in E[y] with gcd(fci, k2) = 1. 



The algorithm ReduceCert in [10] computes a rational 
function Si s.t. 



S = D y (Si) + SiK + 



bk 2 



(4) 



where a, 6 G E[y] satisfy the following conditions: b is the 
squarefree part of the denominator of S, and gcd(o, k2)=l. 
Note that a is not necessarily coprime with bk 2 . As the al- 
gorithm ReduceCert only reduces the shell S, it is referred 
to as the shell reduction. It follows from (4) that 



H = D y (SiT) + —T. 

bk 2 



(5) 



By Theorem 4 in [10], a/6 belongs to E[y] if is integrable. 

Example 4. iet H be the same hyperexponential function 
as in Example 2. Then D y (H)/H has kernel K = y/(y 2 + 1) 
and shell S = l/(y — l) 2 . The shell reduction yields 

S = D y (Si) + SiK + 



(w-i)fe' 

where Si = 1/(1 — y) ana! = J/ 2 + l. TTien i/ can fee decom- 
posed into H = D y (S 1 T)+yT/((y-l)k 2 ), where T=y?y 2 +1. 
By Theorem 4 in [10], H is not integrable. 

On the other hand, it is possible that a in (5) is nonzero 
but H is integrable. 

Example 5. Let H = yexp(y) whose logarithmic derivative 
has kernel 1 and shell y, that is, H = yexp (J ldy) . But H 
is integrable as it is equal to D y (yexp(y) — exp(j/)). 

The shell reduction cannot be directly used to decide hy- 
perexponential integrability. To amend this, the solution 
proposed in [10, Algorithm ReduceHyperexp] was to find 
the polynomial solutions of an auxiliary first-order linear 
differential equation. In the following section, we show how 
this can be avoided and improved. 

4. HERMITE REDUCTION FOR 

HYPEREXPONENTIAL FUNCTIONS 

After the shell reduction described in (5), the denomina- 
tors of shells have been reduced to squarefree polynomials. 
In the rational case, i.e., when the kernel K is zero, the poly- 
nomial a in (4) can be chosen s.t. deg(a) < deg(b), because 
all polynomials are rational integrable. But a hyperexpo- 
nential function with a polynomial shell is not necessarily 
integrable. For example, H = exp (y 2 ) . 

We present a differential variant of [2, Theorem 7] to 
bound the degree of a in (4). The variant leads not only 
to a canonical additive decomposition of hyperexponential 
functions, but also a direct way to decide their integrability. 

4.1 Polynomial reduction 

With Convention 3, we define 

M K = {k 2 D s (p) + ki P \pe E[y]}. 

It is an E-linear subspace in E[j/]. We call Mk the subspace 
for polynomial reduction w.r.t. K. Moreover, define an E- 
linear map (f>K from E[y] to Mk that, for every p G E[y], 
maps p to k 2 D y (p) + kip. We call (f>K the map for polynomial 
reduction w.r.t. K. 

Concerning the subspace Mk and the map 4>k, we have 



Lemma 6. (i) If k 2 D y (g) + kig G E[y] for some g G E(y), 
then g 6 E[y]. (ii) The map cj>K is bijective. 

Proof. Assume that g has a pole. Without loss of generality, 
we assume that the pole is y — and has order m, because 
the following argument is also applicable over the algebraic 
closure of E. Expanding g around the origin yields 
T 

g = h terms of higher orders in y, 

y m 

where r G E\{0}. It follows from k 2 D y (g) + kig G E[y] that 
y = is a pole of 

mr \ ( r 
— + higher terms I + K I V higher terms 

ym + l J \ y 7n 

with order no more than that of K. This implies that y—0 
is a simple pole of K with residue m, which is incompatible 
with K being differential-reduced. The first assertion holds. 

The map <j>K is surjective by its definition. If 0x(p)=O for 
some nonzero polynomial p£E[y], then K equals —D y (p)/p, 
which is nonzero since K ^ 0. So K is not differential- 
reduced, a contradiction. The second assertion holds. ■ 



An E-basis of Mk is called an echelon basis if distinct 
elements in the basis have distinct degrees. Echelon bases 
always exist and their degrees form a unique subset of N. 
Let B be an echelon basis of Mk- Define 

Mk = span E \x l \ i g N and t / deg(/) for all f € b\ . 

Then E[y] = Mk © Mk- We call A4 the standard com- 
plement of Mk- Using an echelon basis of Mk, one can 
reduce a polynomial p to a unique polynomial p G Nk 
s.t. p — p G Mk ■ 

In order to find an echelon basis of M k , we set di — deg ki , 
d 2 =degfc 2 , tk~ - lc(fc 1 )/lc(fc 2 ), and B={<t> K {y n ) \n G N}. 
By Lemma 6 (ii), B is an E-basis of Mk- Let p be a nonzero 
polynomial in E[y]. We make the following case distinction. 

Case 1. di > d 2 . Then 

<Pk(p) = lc(fci) lc(p)j/ dl+dcsp + lower terms. 

So B is an echelon basis, in which deg <f>K (y n ) =di+n for 
all n G N. Accordingly, Mk is spanned by l,y, . . . ,y dl ^ 1 . 

Case 2. di = d 2 — 1 and tk is not a positive integer. Then 

4>k{p)~ (deg(p) lc(fc2)+ lc(fci)) lc(p)j/ dl+dcgp +lower terms. 

(6) 

Since tk is not a positive integer, deg</>x (y n ) = di + n. 
Thus, Mk and Mk have the same bases as in Case 1. 
Case 3. di < d 2 — 1. If deg(p) > 0, then 

4>k(p) = deg(p) lc(fc 2 ) lc(p)y' i2+dcg(p) ~ 1 + lower terms. 

Otherwise, degp = and 4>k(p) = kip. Therefore, B is 
again an echelon basis, in which 

degt/iif (1) = di and deg(j>K(y n ) = d 2 + n — 1 for all n > 1. 



Accordingly, Mk has a basis 1, . . . , y dl ,y , ■ ■ ■ ,y~ 
Case 4- di = d 2 — 1 and tk is a positive integer. It follows 
from (6) that deg0(?/ n ) = di + n if n ^ tk- Furthermore, 
for every polynomial p of degree tk, deg((f>K(p)) is of degree 
less than di + tk- So any echelon basis of Mk does not 
contain a polynomial of degree di + tk ■ Set 

2?' = {>(j/ n )|nGN,n/T A -}. 



d 2 -i 



Reducing <j>(y TK ) by the polynomials in B' , we obtain a 
polynomial r of degree less than di. Note that r is nonzero, 
because B is an E-linearly independent set. Hence, B' U {r} 
is an echelon basis of Mk- Consequently, Mk has an E- 

,y dl -\ y dl+TK }. 



It follows that 



basis |l, y, . . . , y d ^(r)-l t y deg(r)+l 



Example 7. Let K = — 6y 3 /(y 4 + 1), which is differential- 
reduced. Then tk = 6. According to Case 4, Mk has an 
echelon basis {y}u{(n - 6)y n+3 +ny n ~ 1 jnGN, . More- 

over, Mk has a basis {l,y 2 ,y 9 } . 

One can reduce the degree and number of terms of a poly- 
nomial using a subspace of polynomial reduction. 

Lemma 8. With Convention 3, we further let di — dcgfci, 
di — degki, and tk = — lc(fci)/lc(fc2). Let Mk be the sub- 
space for polynomial reduction, and Mk its standard com- 
plement w.r.t. K. Finally, let p be a polynomial in E[j/]. 

(i) If di > di or di = di — 1 and tk ^ Z + , then there 
exists q G Mk s.t. p = q mod Mk and degq < di. 

(ii) If di < di — 1, then there exists q G Mk s.t. p = q 
mod Mk, degq < di and the coefficient of y dl in q 
is equal to zero. 

(Hi) If di=di — 1 and tk£% + , then there exists r G E[y] of 
degree less than di s.t. 



p = sy 



+ r mod Mk for some s G E. 



Moreover, sy dl+TK + r belongs to Mk, and r has at 
most di — 1 terms. 

Proof. The lemma is immediate from the E-bases of Mk 
constructed in the above case distinction. ■ 

The next corollary enables us to derive an order bound on 
telescopers for hyperexponential functions. 

Corollary 9. With the notation introduced in Lemma 8, 
there exists V C {y"\n G N} with \V\ < max(di,d 2 - 1) s.t. 
every polynomial in E[y] can be reduced modulo Mk to an 
E-linear combination of the elements in V '. 

Proof. By the above case distinction, the dimension of Mk 
over E is at most max(di , d% — 1). The corollary follows. ■ 

4.2 Hyperexponential integrability 

With Convention 3, we further assume that the polynomi- 
als a and b are obtained by the shell reduction in (5). So the 
decomposition (5) holds for the present notation. Moreover, 
let Mk be the subspace of polynomial reduction w.r.t. K, 
and Mk its standard complement. 

We are going to determine necessary and sufficient condi- 
tions on hyperexponential integrability. Since gcd(t>, ki)=l, 



a 

bki 



q r 
b ki 



(7) 



where p, q,r G E[y], deg(q) < deg(fe), and deg(r) < deg(ki). 
Using an echelon basis of Mk , we compute u in Mk and v 
in Mk s.t. kip + r = u + v. By the definition of Mk, there 
exists w in E[j/] s.t. u — kiD y (w) + k\w. By (7), we get 



a _ q kiD y (w) + kiw + v 
bki b ki 



Dy(w) + Kw + Y + T~- 
b k 2 



Wi T = D ^ wT)+ {l + Ti^ T - 



(8) 



The process for obtaining (8) is referred to as the polynomial 
reduction for (a/(bki))T w.r.t. K, as it makes essential use 
of the subspaces Mk and Mk- By (8) and (5), 

H = D v ((S 1 +w)T)+(± + ^T, (9) 

which motivates us to introduce the notion of residual forms. 

Definition 10. With Convention 3, we further let f be a 
rational function in E(?/) . Another rational function rGE(y) 
is said to be a residual form of f w.r.t. K if there exist g 
in E(t/) and q, b,p in E[y] s.t. 



f — D y (g) + Kg + r and 



q v 
b + fe' 



where b is squarefree, gcd(b,ki) = 1, degq < degb, and v 
is in the standard complement Mk of the subspace of poly- 
nomial reduction w.r.t. K. For brevity, we say that r is a 
residual form w.r.t K if f is clear from context. 

Residual forms are closely related to the integrability of 
hyperexponential functions. 

Lemma 11. With Convention 3, we further assume that r 
is a nonzero residual form w.r.t. K. Then the hyperexpo- 
nential function rT is not mtegrable. 

Proof. Let Mk be the subspace for polynomial reduction, 
and Mk its standard complement w.r.t. K. By the defini- 
tion of residual forms, there exist b,q G E[j/] with b being 
squarefree and v G Mk s.t. 

q v 

degb > degq, gcd(6, k 2 ) = 1, and r = - + — . (10) 

o ki 

Thus, r can be rewritten as a/ (bki) for some a G E[j/]. Note 
that a is not necessarily coprime with bki. It follows that 



rT 



■ exp 



ki — D y (k 2 
k 2 



dy 



Since (fci — D y (ki)) /ki is differential-reduced and ki,b are 
coprime, (a/b,(k\ — D y (ki))/ki) is indecomposable accord- 
ing to Definition 2 in [10]. By Theorem 4 in [10], a/b is 
in E[j/]. So the denominator of r divides ki, which, together 
with (10), implies that q — 0. Consequently, (v/ki)T is in- 
tegrable. By (2), v = kiD y (s) + fcis for some s G E(y). 
Since v G E[y], v G Mk by Lemma 6 (i). Thus, v = be- 
cause v G Mk- We have that r = 0, a contradiction to the 
assumption that r / 0. ■ 

The existence and uniqueness of residual forms are de- 
scribed below. 

Lemma 12. With Convention 3, we have that the shell S 
has a residual form w.r.t. the kernel K. If a rational function 
has two residual forms w.r.t. K, then they are equal. 

Proof. By (9), S = D y (Si + w) + (Si + w)K + q/b + v/k 2 . 
So q/b + v/ki is a required form. 

Let r = q/b + v/ki and r' = q' /&' + v' /ki be two residual 
forms of a rational function w.r.t. K, where b,b' ,q,q' ,v,v' 
are inE[y], b and b' are squarefree, gcd(6, ki)= gcd(fe', ki) = l, 



degg < degb, degg' < degfo' and t),v' £ Nk- By the defini- 
tion of residual forms, D y (f) + fK + r = D v (f) + f'K + r' 
for some /, /' £ E(i/). It follows that 

D y (,f-f') + {f-f')K + r-r' = 0. 

Hence, (r' — r)T is integrable by (2). Since r — r' is also a 
residual form w.r.t. K, r = r' by Lemma 11. ■ 

Below is the main result of the present section. 

Theorem 13. Let H be a hyperexponential function whose 
logarithmic derivative has kernel K and shell S. Then there 
is an algorithm for computing a rational function h in E(j/) 
and a unique residual form r w.r.t. K s.t. 

H — D y ^/iexp (J K dyj) + rexp (J K dy^j . (11) 

Moreover, H is integrable if and only if r — 0. 

Proof. Let T = exp (J Kdy). Applying the shell reduc- 
tion to H w.r.t. K, we can find a rational function Si, 
and two polynomials a, b s.t. (4) holds. Then we apply 
the polynomial reduction to a/(bk 2 )T to get the residual 
form r = q/b + v/k 2 s.t. (11) holds. 

Suppose that there exists another decomposition 



H — D y (h'T) +r'T 



(12) 



for some h' £ E(t/) and r' is a residual form w.r.t. K. Then 
both r and r' are residual forms of S by (11), (12) and the 
fact H = ST. So r = r' by Lemma 12. 

If r = 0, then H is obviously integrable. Conversely, 
assume that H is integrable. Then rT is also integrable 
by (11). So r = by Lemma 11. ■ 

The reduction algorithm described in the proof of Theo- 
rem 13 has three interesting features. First, it enables us to 
decide hyperexponential integrability immediately. Second, 
it decomposes a hyperexponential function into a sum of an 
integrable one and a non-integrable one in a canonical way. 
Third, it does not need to compute a polynomial solution 
of any first-order linear differential equation. The method 
will be referred to as Hermite reduction for hyperexponen- 
tial functions in the sequel, because it extends all important 
conclusions obtained by Hermite reduction for rational func- 
tions to hyperexponential ones. 

Example 14. Let H be the same hyperexponential function 
as in Example 2. Then K — y/(y 2 + 1) and S = l/(y — l) 2 . 
Set T — \J y? + 1. By the shell reduction in Example 4, 

H = D «{jk T ) + ik T > 

where b — y — 1 and k 2 = y 2 + 1. The polynomial reduc- 
tion yields (y/(bk 2 ))T = D y (-T/2) + (1/(26) + l/(2fe))T. 
Combining the above equations, we decompose H as 



H = D V 



-(, + 1) Wl_ J_ 
2(y -1) J \2b 2k 2 



T. 



Example 15. Consider H — yexp(y) as given in Exam- 
ple 5. Since its logarithmic derivative has kernel K = 1, the 
subspace Mk of polynomial reduction is equal to E[y]. Thus, 
y £ Mk and H is integrable. More generally, Mk ~ E[j/] 
corresponds to the wellknown fact that p(y) exp(y) is inte- 
grable for all p £ E[y] \ {0}. 



5. KERNEL REDUCTION 

Let K = kijki be a nonzero differential-reduced rational 
function in E(y) with gcd(fci,&2) = 1- We may want to 
reduce a hyperexponential function in the form 



_£ 

k r ' 



exp 



j K dy^j for some p £ E[y] and 



m £ 



One way would be to rewrite the above function as 
fci — mD y {k2) 



pexp 



k-i 



dy 



and proceed by polynomial reduction w.r.t. the new ker- 
nel (ki — mD y (k2)) /fej, which is also differential-reduced. 
However, it will prove to be more convenient in Section 6 to 
reduce the given function w.r.t. the initial kernel K. To this 
end, we introduce another type of reduction, based on the 
ideas in [9, 18]. 

Lemma 16. With Convention 3, we letp £ E[y] and m > 1. 
Then there exist pi,pi £ E[y] s.t. 



= Dy 



Pi 



K~ x k 2 



(13) 



Proof. We proceed by induction on m. If m = 1, then tak- 
ing pi = and p2 — p yields the claimed form. Assume 
that m > 1. We first show that there exist p~i,P2 £ EM s.t. 



% 2 \"<2 

which is equivalent to 

p = pi(ki - (m- l)D y (k 2 )) + (D y (pi) +p 2 )fc 2 . 

Since ki/k 2 is differential-reduced, there exist u, v £ E[j/] 
s.t. p = u(ki — (m — l)D y (k 2 )) + vk 2 by the extended Eu- 
clidean algorithm. So we can take pi=u and p 2 =v—D y (u). 
By the induction hypothesis, there exist pi,p 2 £ E[y] s.t. 



>„m-l » 



_£]_] I Pl K | P£ 
„m-2 J "T" u ™~2 ' 



Setting pi = pife +pi and p 2 = p 2 completes the proof. 
With Convention 3, we have 



T^T^Dy 



_pi T \ + El T 

K' 1 ) k 2 



by Lemma 16. This reduction will be referred to as the 
kernel reduction for (jp/k™)T w.r.t. K. 

6. TELESCOPING VIA REDUCTIONS 

Hermite reduction has been used to construct telescopers 
for bivariate rational functions in [5]. The goal of this sec- 
tion is to develop a reduction-based telescoping method for 
bivariate hyperexponential functions. 

6.1 Creative telescoping for bivariate rational 
functions 

We briefly recall the reduction-based method for rational- 
function telescoping in [5]. 

Let F be a field of characteristic zero and ¥(x,y) be the 
field of rational functions in x and y over F. Let D x and D y 
denote the usual derivations d/dx and d/dy, respectively. 



Let ¥(x){D x ) be the ring of linear differential operators 
over F(x). The ring ¥(x){D x ) is a left Euclidean domain 
and its left ideals are principal. For r £ ¥(x,y), the tele- 
scoping problem is to construct a nonzero linear differen- 
tial operator L(x,D x ) £ ¥(x){D x ) s.t. L(x, D x )(r) = D y (s), 
where s £ ¥(x,y). The operator L is called a telescoper 
for r, and s is the corresponding certificate. The set T of 
all telescopers for a given rational function is a left ideal 
of ¥(x){D x ). Any generator of T is called a minimal tele- 
scoper for the given rational function. 

For any i 6 N, rational Hermite reduction (w.r.t. y) de- 
composes D x (r) into D x (r)=D y (si)-\-ai/b, where Si£F(:r,i/) 
and di, b£¥(x)[y] with deg y (ai)< deg y (6), and 6 is square- 
free over ¥(x). Since deg y (ai) is bounded by deg y (b), the se- 
quence {a 4 } ieN is linearly dependent over ¥(x). Assume that 
there exist eo, . . . , e p £ ¥(x), not all zero, s.t. Y2i=o eiCLi = 0- 
Then L := X^=o e »^a; ls a telescoper for r and ~}2i-o ei 3* ' s 
the corresponding certificate. In fact, L is a minimal tele- 
scoper for r if p is the minimal integer s.t. eo,...,e p are 
linearly dependent over ¥(x). This reasoning yields the up- 
per bound deg H (6) on the order of minimal telescopers. 

6.2 Creative telescoping for bivariate hyper- 
exponential functions 

We now apply the Hermite reduction for univariate hyper- 
exponential functions in Section 4 to compute telescopers for 
bivariate hyperexponential functions. 

A nonzero element H in some differential field extension 
of ¥(x, y) is said to be hyperexponential over ¥(x, y) if its log- 
arithmic derivatives D X (H)/H and D y (H)/H are in ¥(x, y). 

Put f=D x (H)/H and g=D y (H)/H. Then D y (f)=D x (g) 
because D x and D y commute. Therefore, it is legitimate to 
denote H by exp(J / dx + gdy). For two hyperexponential 
functions Hi= exp(J /, dx + gi dy), i — 1,2, we have 



H1H2 — exp 



(J (/i + /a)dx + (ffi+fla) 



dy 



In particular, the product of a rational function r £ ¥(x,y) 
and a hyperexponential function H = exp(J f dx + gdy) is 



r# = exp / (/ + D x {r)/r) dx + (g + D y (r)/r) dy 



The following fact is immediate from [10, Lemma 8]. 

Fact 17. Let f and g be rational functions in ¥(x, y) sat- 
isfying D y (f) — D x (g). Then the denominator of f divides 
that of g in ¥(x) [y] . 

For a hyperexponential function H over ¥(x,y), the tele- 
scoping problem is to construct a linear ordinary differential 
operator L(x,D x ) in ¥(x){D x ) s.t. 

L(x,D x )(H) = D y (G) 

for some hyperexponential function G over ¥(x, y). As in the 
rational case, our idea is to apply the Hermite reduction for 
univariate hyperexponential functions w.r.t. y to the deriva- 
tives D X (H) iteratively, and then find a linear dependency 
among the residual forms over ¥(x). 

Lemma 18. Let H = exp(J f dx + gdy) be a hyperexpo- 
nential function over¥(x,y). Let K be the kernel and S the 
shell of g w.r.t. y. Then, for every i £ N, the i-th deriva- 
tive D X (H) can be decomposed into 



where m £ ¥{x,y), T = exp(/(/ - D x (S)/S)dx + K dy) 
and ri £ ¥{x,y) is a residual form w.r.t. K. Moreover, 
let be the denominator of K, b the squarefree part of the 
denominator of S, and Mk the standard complement of the 
subspace for polynomial reduction w.r.t. K. Then 



(15) 



b k 2 

for some qi £ F(x)[t/] with deg y qi < deg y b and Vi £ Mk- 

Proof. We proceed by induction on i. If i = 0, then the 
assertion holds by Theorem 13. 

Assume that D X (H) can be decomposed into (14) and 
assume that (15) holds. Moreover, let f = f — D X (S)/S. 
Consider the (i + l)-th derivative D x + (h). There exists a 
polynomial a in F(x)[j/] s.t. / = ajk^ by D y (f \ = D X (K) 
and Fact 17. A direct calculation leads to 

D X +1 (H) =D y (D x (u lT )) +(^ + + I 7 

\ OK2 K2 

_ ' -qiDxjb) {a - D x (k 2 ))vi " ? , 



l> 2 



Applying the shell reduction to (—qiD x (b) /b 2 ) T and the 
kernel reduction to ((a — D x {k-i))vi/k2) T w.r.t. y, we get 

where wi,W2,pi and p 2 are in F(x)[y]. We then apply poly- 
nomial reduction to ST w.r.t. K, where 

g _ W2 P2 aqi Dx(qi) , D x (vt) 



bk2 k2 bk2 b 



k> 



which leads to 



S = D y ( W )+ W K+{^±L + - V ^- J . 

where w £ ¥(x,y) and q i+1 /b + v i+ i/k2 is the residual form 
of S w.r.t. K. It follows from a direct calculation that 



D?\H) = D y (u i+1 T) + ^1±L + T , 

where m+\ — D x (ui) + u%f + wi/b + pi/fe +w. m 

The main results in the present section are given below. 

Theorem 19. With the notation introduced in Lemma 18, 
we let L = X^=o e iD y with eo, ■ ■ ■ ,e p £ ¥(x), not all zero. 

(i) L is a telescoper for H if and only i/X/i=i eir ' = 0- 

(ii) The order of a minimal telescoper for H is no more 
than deg y (6) + max(deg y (fei), deg y (k 2 ) - 1). 

Proof. We set E = ¥(x) and view that hyperexponential 
functions involved in the proof are over ¥(y). Moreover, 

let u = X)i=o eiUi antl r = ~£2i=o e ' ri - B y ( 14 )' we nave 



L(H) = D y (uT) + rT. 



(16) 



Dl(H) = D v ( Ui T) + nT, 



(14) 



If r = 0, then L is a telescoper by (16). Conversely, assume 
that L is a telescoper of h. Then rT is integrable w.r.t. y 



by (16). Since r is a residual form, it is equal to zero by 
Lemma 11. The first assertion is proved. 

Set A = max(deg y (fci),deg (fe) — 1). Let the residual 
form n = qi/b + be as defined in (14) and (15). By 

Corollary 9, the Wj's have a common set V of supporting 
monomials with \V\ < A. Moreover, deg y (qi) < deg y (b) 
and gcd(6, fe) = 1- Therefore, the residual forms ro, . . . , r p 
are linearly dependent over F(a;) if p > deg y (b) + A. The 
second assertion holds ■ 

Remark 20. By Theorem 19, the first linear dependency 
among the residual forms ro, r±, T2, . . . gives rise to a mini- 
mal telescoper of H. 

Below is an outline of the reduction based telescoping algo- 
rithm for hyperexponential functions, in which the notation 
is that introduced in Lemma 18 is used. 

Algorithm. HermiteTelescoping: Given a bivariate hyper- 
exponential function H = exp(J f dx + gdy) over ¥(x,y), 
compute a minimal telescoper L and its certificate w.r.t. y. 

1. Find the kernel K and shell S of D y (H)/H w.r.t. y. 
Set 6 to be the squarefree part of the denominator of S. 

2. Decompose H into H — D y (uoT) +roT using the Her- 
mite reduction for hyperexponential functions given in 
Theorem 13. If ro = 0, return (1,uqT). 

3. Set p := deg y {b) + max(deg y (fci), deg y (k 2 ) - 1). 

4. For i from to p do 

4.1. Compute (ui,ri) incrementally s.t. 

Di{H) = D y { Ui T)+nT 

by the shell, kernel and polynomial reductions de- 
scribed in Lemma 18. 

4.2. Find rjj G W(x) s.t. £)*._ »foT 3 - = using the al- 
gorithm in [17]. If there is a nontrivial solution, 

return [J2]=o Vi D l, E 3 =o %«i T ) • 

Example 21. Let H = ^x - 2y exp(x 2 y). Then D X {H)/H 
and D y {H)/H are, respectively, 



1 + 4a; 2 y - 8xy 2 
f = 2(x-2y) aUd 9 



-1 + x 3 - 2x 2 y 
x-2y 



Since g is differential-reduced w.r.t. y, g is the kernel and 1 
is the shell of D y (H)/ H w.r.t. y. By Hermite reduction, 



11 = d„ ( ±11) ■ 4r IL 

x K2 



(17) 



Applying D x to the above equation yields 

where r = (-5a; + 8y + Ax z y - 8x 2 y 2 ) /(2x a k 2 .) . The shell, 
kernel and polynomial reduction given in Lemma 18 yields 



D x (H) = D y { 2x2y „ 3 ■// 



3a: 3 - 6 
2x 3 k 2 



H 



(18) 



Combining (17) and (18), we get L = (6 — 3a; 3 ) + 2xD x 
is a minimal telescoper for H and G = (4y — 3x)H is the 
corresponding certificate. 



Remark 22. The algorithm HermiteTelescoping is directly 
based on the proof of Lemma 18. Yet, there is another idea 
for computing a minimal telescoper of H . Namely, we first 
compute a nonzero operator L\ £ ¥(x){D x ) of minimal order 
s.t. Li(H) = D y (Gi) + {p/k2)T for some hyperexponential 
function Gi and polynomial p. Note that such operators al- 
ways exist, because deg y qi in (15) is less than deg y b. Then 
we apply the algorithm HermiteTelescoping to get a minimal 
telescoper L2 for (p/k2)T . In doing so, any rational function 
with denominator b will not appear when we compute L2 ■ It 
turns out that L2L1 is a minimal telescoper of H . An im- 
plementation on this idea is underway. 

6.3 Comparison with the Apagodu-Zeilberger 
bound 

Assume that 



H = uexp (^-) Y\_p t (x,y) c 

^ ' i — 1 



(19) 



where u, r\ , r2, pi, ■ ■ ■ , p m are nonzero polynomials in F[a;, y] 
and ci , . . . , Cm are distinct indeterminates. Theorem cAZ 
in [4] asserts that the order of minimal telescopers for H is 
bounded by 

m 

« := degy (r 2 ) + max (deg y (n), degj / (r 2 )) +^deg H (pO- 1- 

1=1 

Note that H can be viewed as a hyperexponential function 
over F(ci, . . . , c m )(x, y). We now show that a given above is 
no less than the order bound on minimal telescopers for H 
obtained from Theorem 19 (ii). The kernel and shell of the 
logarithmic derivative D y (H)/H axe 



D y ( Pi ) 



K :=D y [ — ) +J2^ " yKt "' and S := u, 



respectively, because K has no integral residue at any simple 
pole, S is a polynomial in F[x,2/], and D y (H)/H is equal 
to K + D y (S)/S. Let K = ki/k 2 with gcd(fei,fe) = 1. A 
direct calculation leads to 

m 

deg y (ki)< deg y (r 1 ) + deg^ (r 2 ) + ^ deg H (p; ) - 1 , 



and 

m 

deg y (k 2 ) < 2deg y (r 2 ) + ^2deg y (pi). 

1=1 

By Theorem 19, the order of minimal telescopers for H is no 
more than max (deg y (fci), deg y (k 2 ) — l), which is no more 
than a by the above two inequalities. 

Indeed, the order bound in Theorem 19 (ii) may be smaller 
than that in Theorem cAZ. 

Example 23. Let H=q c exp(a/q), where a, q are irreducible 
polynomials in ¥[x,y] with deg y (a) < deg y (q), and c is a 
transcendental constant over F. By Theorem cAZ, a min- 
imal telescoper for H has order no more than 3deg y q— 1. 
On the other hand, the kernel and shell of D y (H)/H are 
equal to (D y (a)q—aD y (q)+cqD y (q)) /q 2 and 1, respectively. 
A minimal telescoper has order no more than 2 deg H q — 1 by 
Theorem 19 (ii). 



In general, Christopher's Theorem states that a hyperex- 
ponential function over F(x, y) can always be written as: 



exp ( — ) Y[pi(x,y) Ci , 



(20) 



where u, v, r\,T2 £ F[as, y], Ci is algebraic over F, and pt 
is in ¥(ci)[x, y], i = 1, . . . ,m. A more explicit description 
on (20) can be found in [7]. So H given in (19) is a spe- 
cial instance for hyperexponential functions. In addition, it 
is easier to compute the kernel and shell than to compute 
the decompositions (19) and (20) when a hyperexponential 
function is given by its logarithmic derivatives. 

7. IMPLEMENTATION AND TIMINGS 

We have produced a preliminary implementation of the 
algorithm HermiteTelescoping in the computer algebra sys- 
tem Maple 16. Our Maple code is available from 

http : / / www4 . ncsu . edu/~schen21/HermiteCT . html 

We now compare the performance of our algorithm to the 
Maple function DEtools[Zeilberger] of the telescoping algo- 
rithm in [3] . The examples for comparison are of the form 



■ exp 



where m £ N, p, q, a, b,u,v € Z[x, y] are irreducible and 
their coefficients are randomly chosen. For simplicity, we 
choose A = deg y (p) = deg y (q), fi = deg y (a) = deg y (b), 
and v — deg y (u) = deg y (v). The runtime comparison (in 
seconds) for different examples is shown in Table 1 , in which 

• ZT: the Maple function DEtools[Zeilberger]. 

• HT: our implementation of HermiteTelescoping. 

• order: the order of the computed minimal telescoper. 

• OOM: Maple runs out of memory. 



(A, 




'A 


m) 


ZT 


HT 


order 


(2, 


0, 


2 


1) 


2.23 


2.43 


5 


(2. 


0, 


2 


2) 


2.21 


2.01 


5 


(3. 


0, 


2 


1) 


8.72 


6.64 


6 


(3. 


0, 


2 


2) 


9.38 


6.56 


6 


(6. 


0, 


1 


1) 


45.35 


24.49 


7 


(6. 


0, 


1 


2) 


43.02 


22.91 


7 


(2. 


2, 


2 


1) 


1405.9 


221.5 


9 


(2, 


2, 


2 


2) 


1398.1 


200.34 


9 


(3. 


0, 


3 


1) 


147.92 


47.23 


8 


(3. 


0, 


3 


2) 


151.20 


44.56 


8 


(3. 


3, 





1) 


207.82 


61.10 


8 


(3. 


3, 





2) 


211.70 


58.63 


8 


(3. 


2, 


1 


1) 


304.61 


62.67 


8 


(3. 


2, 


1 


2) 


331.0 


63.61 


8 


(3. 


1, 


3 


1) 


OOM 


534.87 


10 


(3, 


1, 


3 


2) 


OOM 


522.15 


10 



Table 1: Timings (in sec.) were taken on a Mac OS 
X computer with 4Gb RAM and 3.06 GHz Core 2 Duo 
processor. 
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Remark 24. The orders of the computed minimal tele- 
scopers in our experiments are equal to the predicted order 
bounds in Theorem 19 . 



